10.4 Модели вида ARIMA

Прежде чем перейти к рассмотрению модели ARIMA, познакомимся сначала с двумя другими моделями: скользящего среднего и моделью авторегрессии.

Модель скользящего среднего MA(qq)

Модель скользящего среднего порядка qq или просто MA(qq) предполагает следующую зависимость данных:

yt = μ+εt+θ1εt1+...+θqεtq,y_t\ =\ \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + ... + \theta_q \varepsilon_{t-q},

где yty_t — стационарный ряд со средним μ\mu, а εt\varepsilon_t — гауссовский белый шум, то есть εtN(0,σ2)\varepsilon_t \sim \mathcal{N}(0, \sigma^2) и независимы.

По сути наш ряд yty_t выражается через сумму некоторого фиксированного среднего μ\mu, значения белого шума в текущий момент времени εt\varepsilon_t и не более qq предыдущих значений белого шума, домноженных на некоторые коэффициенты, которые являются параметрами модели.

Рассмотрим некоторые свойства модели MA(qq). Как уже было упомянуто выше, ряд yty_t будет являтьcя стационарным со средним μ\mu. Найдем также Dyt\mathsf{D}y_t. Воспользовавшись свойством независимости для εt\varepsilon_t, можем заключить, что

Dyt=(1+θ12++θq2)σ2.\mathsf{D}y_t = \left(1 + \theta_1^2 + \dots + \theta_q^2\right)\sigma^2.

Посчитаем автоковариационную функцию для ряда yty_t, то есть найдем значение cov(yt,yt+τ)cov(y_t, y_{t + \tau}). Легко понять, что если τ>q\tau > q, то cov(yt,yt+τ)cov(y_t, y_{t+\tau}) = 0, т.к. εt\varepsilon_t независимы. Если же τq\tau \leq q, то тогда

cov(yt,yt+τ)=(θτ+θ1θτ+1++θqτθq)σ2.cov(y_t, y_{t+\tau}) = \left(\theta_{\tau} + \theta_1\theta_{\tau + 1} + \dots + \theta_{q-\tau}\theta_{q}\right)\sigma^2.

Записав более компактно, можем получить:

cov(yt,yt+τ)={σ2j=0qτθjθτ+j,τq;0τ>q;cov(y_t, y_{t+\tau}) = \begin{cases} \sigma^2\sum_{j=0}^{q-\tau}\theta_j\theta_{\tau + j}, & \tau \leq q; \\ 0 & \tau > q; \end{cases}

где θ0=1\theta_0 = 1. Из посчитанных значений для дисперсии и ковариационной функции, можете попробовать получить выражение и для автокорреляционной функции. Ее особенностью будет как раз равенство нулю на лаге, превосходящим qq.

Посмотрим на визуализацию:

ma2

Модель авторегрессии AR(pp)

Модель авторегрессии для временного ряда можно записать следующим образом:

yt = α+φ1yt1+...+φpytp+εt,y_t\ =\ \alpha + \varphi_1 y_{t-1} + ... + \varphi_p y_{t-p} + \varepsilon_t,

где yty_t — стационарный ряд, а εt\varepsilon_t — гауссовский белый шум, то есть εtN(0,σ2)\varepsilon_t \sim \mathcal{N}(0, \sigma^2) и независимы. Отметим, что, вообще говоря, для стационарности нужны некоторые условия на коэффициенты φ1,...,φp\varphi_1, ..., \varphi_p.

По сути наш ряд yty_t выражается через сумму некоторого фиксированного числа α\alpha, значения белого шума в текущий момент времени εt\varepsilon_t и не более pp предыдущих значений этого же ряда, домноженных на некоторые коэффициенты, которые являются параметрами модели.

Другими словами, модель AR(pp) — это модель линейной регрессией\textit{линейной регрессией} для которой

  • Таргет: yty_t — значение ряда в момент времени tt
  • Признаки: yt1,...,ytpy_{t-1}, ..., y_{t-p} — значения ряда в предыдущие моменты времени

Введем LL — оператор сдвига, обладающий следующими свойствами:

  • применение LL к ряду дает предыдущее значение этого же ряда: Lyt=yt1.Ly_t = y_{t-1}.
  • применение LL к белому шуму дает предыдущее значение шума: Lεt=εt1.L\varepsilon_t = \varepsilon_{t-1}.
  • применение LL к константе — это константа: Lc=c.Lc = c.

Оператор LL иногда называют также лаговым оператором. Можно рассматривать функции от оператора сдвига, например, кратное применение оператора LL: L2yt=L(Lyt)=L(yt1)=yt2L^2 y_t = L(L y_t) = L(y_{t-1}) = y_{t-2} или L1yt=yt+1L^{-1} y_t= y_{t+1}. Для записей некоторых моделей временных рядов будет удобно использовать лаговый многочлен:

φ(L)=i=1pφiLi\varphi(L) = \sum_{i=1}^p \varphi_i L^{i}

Обратным к оператору φ(L)\varphi(L) называют оператор φ1(L)\varphi^{-1}(L) такой, что:

φ(L)φ1(L)yt=φ1(L)φ(L)yt=yt\varphi(L)\varphi^{-1}(L)y_t = \varphi^{-1}(L)\varphi(L)y_t=y_t

Так, например, для φ<1\lvert \varphi \rvert < 1 можно заключить, что:

11φL=(1φL)1=i=1φiLi\frac{1}{1 - \varphi L} = \left(1-\varphi L \right)^{-1} = \sum_{i=1}^{\infty}\varphi^{i}L^{i}

Рассмотрим модель AR(pp):

yt = α+φ1yt1+...+φpytp+εty_t\ =\ \alpha + \varphi_1 y_{t-1} + ... + \varphi_p y_{t-p} + \varepsilon_t

С помощью оператора сдвига ее можно представить в следующем виде:

a(L)yt = α+εt,a(L) y_t\ =\ \alpha + \varepsilon_t,

где a(z)=1φ1z...φpzpa(z) = 1 - \varphi_1 z - ... - \varphi_p z^p — характеристический полином.

Сформулируем пару важных утверждений:

  • Любой стационарный (в широком смысле) процесс представим в виде MA()MA(\infty), то есть в виде модели скользящего среднего с неограниченным количеством слагаемых (конечное или бесконечное число). Этот результат так же известен как теорема Волда о декомпозиции временного ряда.
  • Модель AR(p)AR(p) задает стационарный временной ряд \Longleftrightarrow все комплексные корни a(z)=0a(z)=0 лежат вне единичного круга.

Приведем пояснение второго утверждения. В самом деле, пусть z1,...,zpz_1, ..., z_p — все его комплексные корни (их ровно pp с учетом кратности), тогда справедливо представление:

a(z)=(zz1)...(zzp)=z1...zp(1zz1)...(1zzp)a(z) = (z-z_1)...(z-z_p) = z_1...z_p \left(1 - \frac{z}{z_1}\right) ... \left(1 - \frac{z}{z_p}\right)

Тогда при представлении временного ряда в виде

yt = α+εta(L)y_t\ =\ \frac{\alpha + \varepsilon_t}{a(L)}

и дальнейшего его разложения на простые дроби возникнут слагаемые вида

εt1Lzj.\frac{\varepsilon_t}{1 - \frac{L}{z_j}}.

Если при этом zjz_j лежит внутри единичного круга или на его границе, то соответствующий ряд будет расходящимся. На самом деле, случай zj=1z_j=1 мы в дальнейшем учтем.

В качестве примера рассмотрим подробнее модель AR(1)AR(1).
Зависимость имеет вид yt=α+φyt1+εty_t = \alpha + \varphi y_{t - 1} + \varepsilon_t, где εtN(0,σ2)\varepsilon_t \sim \mathcal{N}(0, \sigma^2). Для данного ряда можно выписать следующие свойства:

  • Уравнение 1φz=01 - \varphi z = 0, имеет корень λ=1/φ\lambda = 1/\varphi.
  • Тем самым, AR(1)AR(1) стационарен \Longleftrightarrow φ<1\lvert\varphi\rvert < 1. Кроме того, чем меньше φ\varphi, тем предыдущее значение ряда вносит меньший вклад в текущее значение.
  • Если ряд стационарен, то:
    • Eyt=α1φ\mathsf{E} y_t = \frac{\alpha}{1-\varphi}
    • Dyt=σ21φ2\mathsf{D} y_t = \dfrac{\sigma^2}{1-\varphi^2}
    • cov(yt,yth)=φhσ21φ2cov (y_t, y_{t - h}) = \varphi^h \cdot \dfrac{\sigma^2}{1-\varphi^2}.

Разберем первое равенство, остальные получаются аналогично. Возьмем математическое ожидание в уравнении ряда

Eyt=α+φEyt1+Eεt\mathsf{E} y_t = \alpha + \varphi \mathsf{E} y_{t - 1} +\mathsf{E} \varepsilon_t

Поскольку ряд стационарен, то его математическое ожидание не меняется во времени, а для белого шума математическое ожидание равно нулю. Тем самым мы получаем уравнение на m=Eytm = \mathsf{E} y_t, откуда следует доказываемая формула.

Таким образом, в зависимости от значения φ\varphi мы можем получить следующие результаты:

  • Если φ<1\lvert\varphi\rvert < 1, то yt=μ+j=0φjεtjy_t = \mu + \sum\limits_{j = 0}^{\infty}\varphi^j\varepsilon_{t - j} — представление ряда в виде MA(\infty).
  • Если φ=1\lvert\varphi\rvert = 1, то AR(1)AR(1) — это случайное блуждание.
  • Если φ>1\lvert\varphi\rvert > 1, то AR(1)AR(1) — экспоненциально растущий процесс.

Посмотрим на визуализацию.

ar1

  • В первом случае мы имеем модель yt=0.5yt1+εty_t = - 0.5 y_{t - 1} + \varepsilon_t, отрицательный коэффициент является следствием больших колебаний ряда.

  • Во втором случае модель yt=0.9yt1+εty_t = 0.9 y_{t - 1} + \varepsilon_t, большой положительный коэффициент делает ряд менее шумным.

  • В третьем случае показано несколько рядов вида случайного блуждания yt=yt1+εty_t = y_{t - 1} + \varepsilon_t, что соответствует случаю φ=1\varphi=1.

  • В четвертом случае показан экспоненциальный процесс yt=1.1yt1+εty_t = 1.1 y_{t - 1} + \varepsilon_t, на графике шум уже не заметен из-за масштаба.

На немного вернемся к модели MA(qq). Чуть выше мы выяснили, что при некоторых условиях на коэффиценты φ\varphi временной ряд модели AR(pp) будет стационарным, а значит имеет представление в виде MA(\infty). На самом деле, модель скользящего среднего порядка qq тоже можно представить с помощью оператора LL следующим образом:

yt = μ+εt+θ1εt1+...+θqεtq = μ+b(L)εty_t\ =\ \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + ... + \theta_q \varepsilon_{t-q}\ =\ \mu + b(L)\varepsilon_t

где b(z)=1+θ1z+...+θqzqb(z) = 1 + \theta_1 z + ... + \theta_q z^qхарактеристический многочлен. Для простоты изложения пусть μ=0\mu = 0. Важным при такой записи оказывается понятие обратимости, то есть представления в виде

εt=b1(L)yt,\varepsilon_t = b^{-1}(L)y_t,

которое означает, что ряд можно представить в виде бесконечной авторегрессионной модели.
Здесь, как и в рассуждениях выше, можно заключить, что временной ряд yty_t обратим, если все комплексные корни b(z)=0b(z) = 0 лежат вне единичного круга.

Модель ARMA(p,qp, q)

Модель ARMA(p,qp, q) по сути является суммой моделей AR(p)AR(p) и MA(q)MA(q), иначе говоря, модель есть сумма нескольких предыдущих значений ряда и нескольких предыдущих значений белого шума с некоторым коэффициентами.

yt = α+φ1yt1+...+φpytp+εt+θ1εt1+...+θqεtqy_t\ =\ \alpha + \varphi_1 y_{t-1} + ... + \varphi_p y_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + ... + \theta_q \varepsilon_{t-q}

Эквивалентную запись ряда в терминах оператора сдвига можно получить, рассмотрев два многочлена

a(L)yt=α+b(L)εta(L) y_t = \alpha + b(L) \varepsilon_t

или

yt=μ+b(L)a(L)εt,y_t = \mu + \frac{b(L)}{a(L)} \varepsilon_t,

где a(z)=1φ1z...φpzp,a(z) = 1 - \varphi_1 z - ... - \varphi_p z^p, и b(z)=1+θ1z+...+θqzqb(z) = 1 + \theta_1 z + ... + \theta_q z^q.
Заметим, что во втором представлении константа α\alpha заменена на μ=Eyt\mu = \mathsf{E} y_t. На самом деле, стационарность такого ряда будет определяться только его AR(pp) компонентой, то есть значениями коэффициентов φ\varphi, так ряд в модели MA(qq) всегда является стационарным.

Модель ARIMA(p,d,qp, d, q)

Модель ARIMA(p,d,qp, d, q) — это расширение моделей типа ARMA на нестационарные временные ряды, которые однако могут стать стационарным после применениея процедуры дифференцирования ряда. Модель ARIMA(p,d,qp, d, q) для ряда yty_t определяется как модель ARMA(p,qp, q) для ряда разностей порядка dd ряда yty_t.

  • Разность порядка 1: ytyt1=(1L)yty_t - y_{t-1} = (1 - L) y_t.
  • Разность порядка 2: (1L)2yt=(1L)(ytyt1)=(ytyt1)(yt1yt2)=yt2yt1+yt2(1 - L)^2 y_t = (1 - L) (y_t - y_{t-1}) = (y_t - y_{t-1}) - (y_{t-1} - y_{t-2}) = y_t - 2y_{t-1} + y_{t-2}.

Получаем формулу модели ARIMA:

a(L)(1L)dyt=α+b(L)εta(L) (1 - L)^d y_t = \alpha + b(L) \varepsilon_t

или

(1L)dyt=μ+b(L)a(L)εt.(1 - L)^d y_t = \mu + \frac{b(L)}{a(L)} \varepsilon_t.

То есть многочлен a~(z)=a(z)(1z)d\widetilde{a}(z) = a(z) (1-z)^d имеет dd единичных корней.
Тем самым такая модель позволяет учесть нестационарности, в частности, тренд.

В качестве примера рассмотрим процесс случайного блуждания:

yt=yt1+εt,y_t = y_{t-1} + \varepsilon_t,

где εt\varepsilon_t — белый шум. Как уже упомяналось ранее, такой ряд не является стационарным. Однако, если мы применим операцию дифференцирования, то можем перейти к новому, уже стационарном ряду yt=ytyt1y'_t = y_t - y_{t-1}, который можно записать в виде:

yt=εty'_t = \varepsilon_t

Частичная автокорреляция

Для модели скользящего среднего порядка qq мы выяснили, что значения автокорреляционной функции для такого ряда оказывается равной нулю после лага qq. Эта особенность позволяет использовать автокорреляционную функцию для определения порядка модели скользящего среднего. Возникает разумный вопрос, как оценить порядок pp для модели AR(pp)? Здесь оказывается полезным понятие частичной (частной) автокорреляционной функции.

Частичная автокорреляция (PACF) — корреляция ряда с собой после снятия линеной зависимости от промежуточных значений ряда. Иначе говоря, мы хотим как-то учесть опосредованного влияние промежуточных значений ряда и оценить непосредственное влияние ytτy_{t - \tau} на yty_t. Чуть более формально частичную автокорреляцию можно записать следующим образом:

γτ={corr(yt+1,yt),τ=1;corr(yt+τyt+τh1,ytytτ1),τ2,\gamma_{\tau} = \begin{cases} corr(y_{t+1}, y_t), & \tau=1; \\ corr\left(y_{t+\tau} - y_{t+\tau}^{h-1}, y_t - y_t^{\tau-1}\right), & \tau\geqslant2, \end{cases}

где ytτ1y_t^{\tau-1} — линейная регрессия на yt1,yt2,...,yt(τ1)y_{t-1}, y_{t-2}, ..., y_{t - (\tau-1)}:

  • ytτ1=φ1yt1+φ2yt2+...+φτ1yt(τ1)y_t^{\tau-1} = \varphi_1 y_{t-1} + \varphi_2 y_{t-2} + ... + \varphi_{\tau-1} y_{t - (\tau-1)}
  • yt+ττ1=φ1yt+τ1+φ2yt+τ2+...+φτ1yt+1y_{t+\tau}^{\tau-1} = \varphi_1 y_{t+\tau-1} + \varphi_2 y_{t+\tau-2} + ... + \varphi_{\tau-1} y_{t+1}

Пример для τ=2\tau=2:

γ2=corr(yt+2φ1yt+1,ytφ1yt1)\gamma_{2} = corr\left(y_{t+2} - \varphi_1 y_{t+1}, y_t - \varphi_1 y_{t-1}\right)

где φ1\varphi_1 — МНК-оценка в модели yt=φyt1y_t = \varphi y_{t-1}.

Можно показать, что значение частиной автокорреляции для модели авторегресии AR(pp) будет ненулевой для лагов τp\tau \leq p и равняться нулю для лагов τ>p\tau > p. Имеет место быть полная аналогия с автокорреляционной функцией и моделью MA(qq). Таким образом, исследование поведения автокорреляционной и частичной автокорреляционной функции может быть использовано для определения порядка qq модели скользящего среднего и порядка pp модели авторегрессии соответсвтенно.

Оценка коэффициентов в ARIMA

Пусть гиперпараметры p,d,qp, d, q фиксированы.&tab;В предположении, что εt\varepsilon_t — гауссовский белый шум, в нашей модели мы можем выписать функцию правдоподобия Ly(θ,φ,α)=pθ,φ,α(y1,...,yT),L_y(\theta, \varphi, \alpha) = p_{\theta, \varphi, \alpha}(y_1, ..., y_T), где pθ,φ,α(a1,...,aT)p_{\theta, \varphi, \alpha}(a_1, ..., a_T) — соместная плотность. Из-за того, что εt\varepsilon_t имеют нормальное распределение, она будет иметь разумный вид. Соответственно, в качестве оценок параметров берется оценка максимального правдоподобия.

Для поиска начальных приближение для параметров pp и qq воспользуемся автокорреляционной и частичной автокорреляционной функцией.

  • Начальное приближение pp: последний значимый пик у PACF.
  • Начальное приближение qq: последний значимый пик у ACF.

Далее обычно используется поиск по сетке вокруг подобранных значений, минимизируя информационный критерий:

  • AIC=2+2(p+q+1)AIC = -2\ell^* + 2(p+q+1) — критерий Акаике;
  • AICc=2+2(p+q+1)(p+q+2)Tpq2AIC_c = -2\ell^* + \frac{2(p+q+1)(p+q+2)}{T-p-q-2} — критерий Акаике (короткие ряды);
  • BIC=2+(logT2)(p+q+1)BIC = -2\ell^* + (\log T - 2)(p+q+1) — Байесовский информационный критерий или критерий Шварца,

где =lnLy(θ^,φ^,α^)\ell^* = \ln L_y\left(\widehat{\theta}, \widehat{\varphi}, \widehat{\alpha}\right) — логарифм функции правдоподобия, TT — длина временного ряда.

Приведем некоторый план при применению модели ARIMA для прогнозирования временных рядов.

  1. Анализ выбросов: замена нерелевантых выбросов на NA или усреднение по соседним элементам.

  2. Стабилизация дисперсии (преобразования).

  3. Дифференцирование, если ряд не стационарен.

  4. Выбор пилотных pp и qq по PACF и ACF.

  5. Вокруг этих параметров подбираем оптим. модель по AICAIC/AICcAIC_c.

  6. Пошаговое построение прогноза:
    — для tTt \leqslant T: εtε^t=yty^t\varepsilon_t \Longrightarrow \widehat{\varepsilon}_t = y_t - \widehat{y}_t;
    — для t>Tt > T: εt0\varepsilon_t \Longrightarrow 0;
    — для t>Tt > T: yty^ty_t \Longrightarrow \widehat{y}_t.

  7. Построение предсказательного интервала:

    — если остатки модели нормальны и гомоскедастичны (дисперсия постоянна), то строится теоретический предсказательный интервал

    σ^2(h)=σ^2(1+i=1h1ψ^i2)\widehat \sigma^2(h) = \widehat \sigma^2 \left(1 + \sum\limits_{i = 1}^{h - 1} \widehat{\psi}_i^2\right)

    где hh — горизонт прогнозирования, σ^2\widehat\sigma^2 — оценка на дисперсию шума εt\varepsilon_t, ψ^i\widehat{\psi}_i — коэф. для ряда при его представлении в виде бесконечного процесса скользящего среднего. И σ^2\widehat\sigma^2, и ψ^i\widehat{\psi}_i могут быть выражены через оценки на параметры φ\varphi и θ\theta.
    — иначе интервалы строятся с помощью бутстрепа.

Модели SARIMA и ARIMAX

Рассмотрим некоторые расширение модели ARIMA. Обобщение модели ARIMA на ряды с наличием сезонной составляющей назвается SARIMA. Пусть ss — известная сезонность ряда. Добавим в модель ARIMA(p,d,qp, d, q) компоненты, отвечающие за значения в предыдущие сезоны. Тогда модель SARIMA (p,d,q)×(P,D,Q)s(p, d, q)\times (P, D, Q)_s может быть записана следующим образом:

(1L)d(1Ls)Dyt=μ+b(L)B(Ls)a(L)A(Ls)εt,(1 - L)^d (1 - L^s)^D y_t = \mu + \frac{b(L) B(L^s)}{a(L) A(L^s)} \varepsilon_t,

где

a(z)=1φ1z...φpzp,a(z) = 1 - \varphi_1 z - ... - \varphi_p z^p,

b(z)=1+θ1z+...+θqzq,b(z) = 1 + \theta_1 z + ... + \theta_q z^q,

A(z)=1φ1szφPzP,A(z) = 1 - \varphi_1^s z - \dots - \varphi_P z^P,

B(z)=1+θ1sz++θQszQ.B(z) = 1 + \theta^s_1 z + \dots + \theta_Q^s z^Q.

Параметр сезонного дифференцирования DD, а также параметры P,QP, Q подбираются из тех же соображений, что и для p,d,qp, d, q, но только с поправкой, что делается это с учетом сезонности ss. ARIMAX — обобщение модели ARIMA, которая учитывает некоторые экзогенные факторы. Пусть xtRnx_t \in \mathbb{R}^n — ряд регрессоров, известный до начала прогноза.

Простой вариант:

(1L)dyt=μ+i=1nβia(L)xti+b(L)a(L)εt(1 - L)^d y_t = \mu + \sum_{i=1}^n \frac{\beta_i}{a(L)} x_t^i + \frac{b(L)}{a(L)} \varepsilon_t

Общий случай:

(1L)dyt=μ+i=1nui(L)vi(L)xti+b(L)a(L)εt(1 - L)^d y_t = \mu + \sum_{i=1}^n \frac{u_i(L)}{v_i(L)} x_t^i + \frac{b(L)}{a(L)} \varepsilon_t

Пример: xt=I{в момент времени t праздник}x_t = I \{\text{в момент времени t праздник}\}

Вышеуказанные модели можно объединить и получить SARIMAX (p,d,q)times(P,D,Q)s(p, d, q)\\times (P, D, Q)_s:

(1L)d(1Ls)Dyt=μ+i=1nui(L)vi(L)xti+b(L)B(Ls)a(L)A(Ls)εt(1 - L)^d (1 - L^s)^D y_t = \mu + \sum_{i=1}^n \frac{u_i(L)}{v_i(L)} x_t^i + \frac{b(L) B(L^s)}{a(L) A(L^s)} \varepsilon_t

Проведем аналогию с линейной регрессией. Это линейная по признакам модель, в которой

  • Отклик: yty_t — значение ряда в моменты времени tt
  • Признаки:
    • yt1,...,ytpy_{t-1}, ..., y_{t-p} — значения ряда в предыдущие моменты времени
    • Значение ряда за предыдущие сезоны
    • Значения признаков в предыдущие моменты времени
    • Значения признаков в предыдущие сезоны
  • Ошибка: сумма шума за предыдущие моменты времени и предыдущие сезоны.
Чтобы добавить в заметки выделенный текст, нажмите Command + E

Пройдите квиз по параграфу

Чтобы закрепить пройденный материал
Предыдущий параграф10.3. Аналитика временных рядов
Следующий параграф10.5. Задача ранжирования

Вступайте в сообщество хендбука

Здесь можно найти единомышленников, экспертов и просто интересных собеседников. А ещё — получить помощь или поделиться знаниями.